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We discuss generic limits posed by the trap in atomic systems on the accurate determination of 
critical parameters for second-order phase transitions, from which we deduce optimal protocols to 
extract them. We show that under current experimental conditions the in-situ density profiles are 
barely suitable for an accurate study of critical points in the strongly correlated regime. Contrary 
to recent claims, the proper analysis of time-of-fight images yields critical parameters accurately. 
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Recent breakthroughs in cold atom lattice experiments 
have paved the road to reliable studies of condensed 
matter physics by emulating its intractable models. Al- 
though the identification of phases in present experiments 
is nearly universally accepted, the mesoscopic size and 
the presence of a curved potential have long plagued our 
understanding of criticality in optical lattice experiments. 
This has sparked vivid debates [U-IH for even the sim- 
ple normal-to-superfluid (N-SF) U(l) symmetry break- 
ing transition for interacting bosons. In order to be reli- 
able quantum emulators, cold atom experiments need a 
generic protocol for studying critical behavior, and they 
need to do so without invoking numerical simulations. 

The original experiments in optical lattices @ were 
based on time-of-flight (TOF) images and the emergence 
of a 'sharp' peak in the momentum distribution. The ap- 
proach is adequate for identifying superfluid and normal 
phases, but it lacks an explicit procedure for locating the 
phase boundary, since peaks already develop in the nor- 
mal phase due to an increasin g co rrelation length 
It has been proposed in Ref. [ll| that one should trace 
the evolution of the peak shape, and use the increase in 
the peak width as the transition signature. The peak 
width keeps decreasing across the transition point how- 
ever, and this leaves room for subjective data analysis 
and thus increased error bars. Monitoring both the peak 
width and the number of particles in the low-momentum 
peak has been used recently in measuring the N-SF tran- 
sition temperature and its suppression on approach to 
the Mott insulating phase 

Difficulties with the analysis of TOF images and ad- 
vances in single-site resolution detection methods 12i.ll3| 
prompted the authors of Ref. 0| to propose that finite- 
temperature critical points of strongly correlated quan- 
tum models emulated by an optical lattice experiment 
can be generically deduced from cusps in the derivative 
of the density profile of atoms in the trap with respect to 
the external potential, n(r) = —dn(r)/dV(r). Within the 
local density approximation (LDA) approach this quan- 
tity coincides with the compressibility dn/d^i. However, 
the proposal of Ref. [3| was shown not to work under 



realistic experimental conditions [14 1. 

The goal of this Letter is to answer qualitatively and 
quantitatively whether determining critical parameters 
in a trapped system using available experimental tech- 
niques is feasible. First, we employ the theory of finite- 
size effects in critical phenomena discussed recently for 
trapped systems by Campostrini and Vicari [15[ to set 
theoretical limits on the accurate location of phase tran- 
sition points in a given trap for a generic, strongly cor- 
related system. Second, we demonstrate that LDA vio- 
lations in density profiles indicate critical behavior, not 
LDA itself. The calculations show, however, that experi- 
mental data should have extremely small error bars, well 
below what is currently available. Taking the numerical 
derivative of the density profiles is of no help, and is al- 
ways either too noisy or lacking features associated with 
criticality. Third, we demonstrate that detection meth- 
ods coupling directly to critical modes, e.g., TOF images 
in the case of the N-SF transition, do allow one to reach 
the theoretical limit. From the analysis of the central 
peak shape width we construct a quantity with a sharp 
minimum at the critical point. Our considerations apply 
to any scale-invariant second-order phase transition. 

We simulate the 3D Bose-Hubbard model, 
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where t is the matrix element describing the hopping of 
bosons between nearest-neighbor sites, U is the on-site 
repulsion, Vi is the trapping potential, and m = b\b i is 
the on-site occupation number in terms of bosonic cre- 
ation and annihilation operators. In what follows, we use 
t as the unit of energy. 

Let the phase transition occur around the chemical po- 
tential /i(r c ) = /Lt c , where \x c is the critical point in a ho- 
mogeneous system, and r c is situated away from the trap 
center. By the very nature of second-order phase transi- 
tions, characterized by a divergent correlation radius £ at 
fi c , there exists a finite shell centered at r = r c in which 
LDA fails. Indeed, LDA implies quasi-homogeneity of 
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the system when the change of thermodynamic proper- 
ties of the system is negligible at the distance of the order 
of correlation radius. Clearly, the quasi-homogeneity is 
preserved as long as |£(r ± £(r)) — £(r)| <C £(r), and is 
violated otherwise, the latter inevitably happening when 
\ r ~ r c\ where we have used the shorthand nota- 

tion = £(M( r ))- The condition 



,1 ~ £(r.) = 



(2) 



defines the position of the boundaries r* of the LDA vio- 
lated region, and the largest possible correlation radius. 
The value of controls the rounding effects, similar to 
a finite-size system, leading to an intrinsic uncertainty in 
determining critical parameters: Sfj. c oc where v is 

the correlation length exponent. To cast all relations in 
dimensionlcss form we introduce a typical scale for the 
chemical potential, fio, which we choose such that when 
jj, is changed by (io the system density changes by a fac- 
tor of two. We also introduce the healing length d as the 
typical length scale. Then, ((fytc/Mo) ~ (d/t;*) 1 ^- 
In a constant potential gradient, Eq. ([2j implies 
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d ~ \dVfi J 
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Assuming spherically symmetric harmonic confinement, 
we have V/Lt = mui 2 r c . By introducing the typical cloud 
size using muS 2 R 2 = /io, the final estimate for Sfi c can be 
identically rewritten as 
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Obviously, the best accuracy is found when the critical 
point is located at the trap center where the chemical po- 
tential gradient is zero. In this case the derivation has to 
be repeated along the same lines using 5jjL c ~ moj 2 £ 2 for 
the change of the chemical potential over the correlation 
length. The self-consistent solution then takes the form 



R 



2v/(2+v) 



Mo 



R 



2/(2i/+l) 



(5) 



reproducing the result of Ref. [15j for scaling of the corre- 
lation length with the trap size. For the N-SF transition 
(y = 0.6715) in the strongly correlated Bose lattice sys- 
tem, the size is R/d ~ 100, which leads to the conclusion 
that the critical point may, in principle, be determined 
with a relative accuracy of a few percent. Note that this 
estimate is assuming that system properties are known 
exactly, i.e., it does not take into account experimental 
noise which we consider next. 

Density profiles n(r) can be measured with single- 
site resolution [13, 0], and, under ideal experimen- 
tal conditions, further averaged over about a hundred 
shots. Unfortunately, the density and its derivative are 



continuous functions across the N-SF transition point. 
More precisely, the critical contribution to compressibil- 
ity, k = dn/dfi, is of the form A(\[i — yU c |//j, ) Q , where 
a = dv — 2 = 0.015 with slightly different amplitudes 
A on the two sides of the transition. For the U/t = 24, 
T/t = 2.4 example considered below the value of A is 
about 0.5/t. To understand the effect of critical fluctua- 
tions on density profiles we invert the definition of com- 
pressibility and calculate the deviation of the density at 
point r c from the thermodynamic limit value n c as 



n(r c ) 
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dx 



Mo 



(6) 



This equation accounts for LDA violations due to the 
difference between the thermodynamic expression for the 
critical part of compressibility and its maximal value in 
a system of size . Given the small value of a for most 
continuous phase transitions we can write 



n(r c ) = n c + 
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which allows us to estimate the size of the effect in realis- 
tic traps using Eqs. ([4]) and {5j. The largest LDA viola- 
tion is expected at the trap perimeter where the chemical 
potential gradients are large. 

Since the critical value of the density is not known 
a priori, one should measure density profiles with dif- 
ferent global chemical potentials Mi( r = 0)j or in traps 
with different confinement frequencies Wj, convert them 
to n(fi,i) curves, and plot them in the same figure. The 
expectation is that densities measured at different chem- 
ical potential gradients will show the largest deviation 
from each other in the critical region. Unfortunately, the 
size of the effect is supposed to be rather small in the 
strongly correlated regime where Afi is not large, see 
Eq. ([7]). The best possible contrast can be found when 
the profile n(/i, 2) with the critical region in the center is 
subtracted from the profile with the critical region in the 
perimeter n(/x, 1), yielding Sn(fi) = n{fi, 1) — n{fi, 2), and 
thus using the system with the smallest size of the criti- 
cal region as the reference curve. From Eqs. (J4j) , (JSJ) , and 
d?]) we deduce that Sn(fi c ) obtained for profiles with the 
critical point at the trap perimeter and the trap center 
is about 



5n{nc) ~ aAfi 



d ^ 1/(1+1/) ( d x 2/(2^+1) 
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which for the N-SF transition with size R/d = 100 is as 
small as 4 x 10~ 4 A/io, i.e., it is about a few tenths of a 
percent only! In Fij 
Carlo simulation 111 



i.e., 

. 1] we show the result of the Monte 
of a trapped system for U/t = 



24, T/t = 2.4, and mu 2 a 2 /2t = 0.001, where a is the 
lattice constant. As expected, the largest LDA violations 
are observed for the curve with the critical point closest 
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to the trap center. From the position of the maximum 
we deduce that fx c = —2.3 ± 0.4, and n c = 0.37 ± 0.05. 

It is thus possible, in principle, to determine critical 
points with an accuracy of about 10-20% by examining 
LDA violations in density profiles provided they are mea- 
sured with extremely high accuracy at the trap center. 
[Note that we implicitly assume that experiments are able 
to determine the system temperature and chemical po- 
tential by other methods, e.g., using the equation of state 
at the trap perimeter, entropy matching, etc.] It is evi- 
dent from Fig.[T]that density profiles have to be measured 
with a relative accuracy of 0.001. As far as we know such 
accuracy is not available in current experiments. The 
noise in the density at the trap center is determined by 
the number of experimental runs because in one mea- 
surement Srii ~ m ~ 1. Optimistically, the error bars on 
radial density are about er(0) <~ 0.1 at the trap center, 
quickly reducing to a(r) ~ a(0)d/4r at a distance r due 
to radial averaging. Hypothetically, it would be possible 
to have data in three-dimensional systems with accuracy 
of 0.001 at a distance r/d > 25 with signal-to-noise ra- 
tion close to unity if all other experimental uncertainties 
are eliminated. 
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FIG. 2: (Color online). Local compressibility across the crit- 
ical point at jj,c = —2.037 (marked with an arrow). In the 
upper panel the critical region is close to the trap center. In 
the lower panel the critical region is close to the trap perimeter 
and is subject to a larger chemical potential gradient. Com- 
pressibility data were deduced from curves shown in Fig. [T] 
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FIG. 3: (Color online). Momentum distributions and the 
analysis of the low-momentum peak shape as described in the 
text. Here, s = 2. 



FIG. 1: (Color online). Density profiles in the vicinity of the 
critical point (left panel) demonstrating violations of the LDA 
approximation in the presence of finite chemical potential gra- 
dient. In the right panel we plot the difference between the 
two curves. The critical point fi c = —2.037 is marked by an 
arrow. 

Our final note concerns density profile derivatives k ~ 
—dn/dV(r). It is clear from the beginning that measures 
based on integrals of k, i.e., the LDA violations discussed 
above, contain all the information and are more sensitive 
because they efficiently average out the noise in numer- 
ical derivatives. Indeed, in the upper panel of Fig. [2] we 
show that n(r) deduced from the profile giving the largest 
signal in Fig. [1] is meaningless. Derivatives taken away 
from the trap center are less noisy, but they also lack 
any signature of the critical point due to the rounding 
of critical singularities by large chemical potential gradi- 
ents, see the lower panel in Fig. [2] In general, LDA viola- 



tions are a direct consequence of a diverging length scale 
while rounded peaks in derivatives are not necessarily 
originating from criticality and may thus be misleading, 
especially when searching for novel phases. 

There is thus no viable alternative to directly studying 
critical modes. For the N-SF transition they are encoded 
in the low-momentum part of the distribution n(k) mea- 
sured in TOF images. Here, we discuss the best case for 
the accuracy when the phase transition starts in the mid- 
dle of the trap. In a homogeneous system, on approach 
to the critical point from the normal phase the peak in 
the momentum distribution around k = is described by 



i(k) 



£2-77 
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for 



£ _1 < k < or 1 
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where rj is the correlation function exponent (for most 
critical points 77 is positive and small; for the N-SF tran- 
sition rj = 0.038). In a trap, we have to account for 
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the fact that £ depends on the position in space since 
for a given detuning at the trap center, A/x(T), we have 
£(r) ~ [A/i(T) + mw 2 r 2 /2]" 1 ' at some distance r. Al- 
though shells away from the center have smaller £ (in 
the normal phase) they have an increasingly larger vol- 
ume. Curiously, in a wide trap the peak amplitude, 
P = n(k = 0), at the critical point is determined by 
regions with short correlation length because the LDA 
integral for P 
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diverges at large r meaning that the major contribution 
comes from a region with £ ~ d, a posteriori justifying 
the use of LDA for evaluating the peak amplitude. 

Thus, at T > T c the TOF image remains relatively 
broad and the critical region at the trap center is not 
dominating in the signal amplitude in the theoretical 
limit of a very broad trap. [ As we show in Fig. [3l it 
remains a sizable 50% of the total signal for a realistic 
trap with R/d ~ 100 ]. Past the critical point, the low- 
momentum part of the distribution develops the conden- 
sate peak smeared out by the presence of the trapping 
potential and the finite size of the SF domain superim- 
posed on the broader peak representing the critical shell 
separating normal and supcrfiuid phases. 

The all-decisive question is how to "read out" the criti- 
cal point. Our proposal is based on monitoring the devel- 
opment of the bimodal structure at low momenta across 
the transition point. We note that the shape (in contrast 
to the value of the signal) of n(k) at small momenta is 
determined by the critical region at the center because 
the first (and also higher-order) derivative, dn(k)/dk, is 
described by an integral over d 3 r which diverges at the 
lower limit! This suggests the following protocol for an- 
alyzing the data: 

(i) find the momentum k max (T) at which the absolute 
value of the first derivative dn(k)/dk has a maximum; 

(ii) construct the "amplitude of the critical signal" as the 
difference P C {T) = n(k = 0) — n(k max ); 

(iii) plot Q(T) = P c (T)fc^ ax (T), with some exponent 
s > 2 — 77, as a function of temperature to "read out" 
the critical point from the minimum in Q(T). 

By construction, P c oc ^ 2 ~ n and fc max oc l/£, and thus 
the choice of s is ensuring that Q(T) is decreasing on the 
normal side of the transition as T approaches T c . On 
the superfiuid side, the formation of the condensate at 
the trap center results in a sharp increase of the peak 
amplitude n(k = 0) and thus Q(T). The critical point 
is then estimated from the position of the minimum in 
Q(T). In Fig. [3] we present results of Monte Carlo simula- 
tions for our reference system, U/t = 24, done at various 
temperatures in a trap with fj,(r = 0) = —2.037, and 
mw 2 a 2 /2 = 0.0033i. The original data for n{k) shown in 
the left panel were analyzed to produce the Q(T) curve 
shown in the right panel. The position of the minimum 



at T/t = 2.4 gives the correct value of the critical tem- 
perature with an accuracy of ~ 5% limited by the tem- 
perature grid; otherwise the error bars are obtained from 
the FWHM of 1/Q. 

If the momentum resolution of TOF images is better 
than the limit 7r/£„, set by Eq. ([5]), including finite time- 
of- flight [l8| , optical broadening and resolution effects [j| , 
then we believe the same accuracy would be possible ex- 
perimentally For columnar images, the exponent s can 
be reduced to > 1 — 77. 

In conclusion, critical points in strongly correlated 
states of trapped atomic systems can be measured with 
relative accuracy of a few percent. The most sensitive 
signal-to-noise probes are based on critical modes with 
diverging correlation functions. 
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